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& I ABSTRACT 

d ■ Millisecond and binary pulsars are the most stable natural frequency standards. They 

can be applied to a number of principal problems of modern astronomy and time- 
c/2 ■ keeping metrology including the search for stochastic gravitational wave background 

in the early universe, testing General Relativity and establishing of new ephemeris 
time scale. The full exploration of pulsar properties requires obtaining proper unbi- 
ased estimates of the spin and orbital parameters, a problem which deserves a special 
investigation. These estimates depend essentially on the random noise component be- 
ing revealed in the residuals of time of arrivals (TOA) and having a different physical 
origin. In the present paper, the influence of low- frequency ("red") timing noise with 
spectral indices from 1 to 6 on TOA residuals, variances, and covariances of estimates 
of measured parameters of single and binary pulsars are studied. In order to determine 
their functional dependence on time, an analytical technique of processing of obser- 
vational data in time domain is developed. Data processing in time domain is more 
informative since it takes into account both stationary component of noise and its 
i non-stationary countepart. Data processing in frequency domain is valid if and only if 

Qs, [ the noise is stationary. Our analysis includes a simplified timing model of a binary pul- 

*"^» i sar in a circular orbit and procedure of estimation of pulsar parameters and residuals 

' under the influence of red noise. We reconfirm, in accordance with results of previous 

" ^ 1 authors, that uncorrelated white noise of errors of measurements of TOA brings on 

gradually decreasing residuals, variances and covariances of all parameters. On the 
other hand, we show that any low-frequency, correlated noise of terrestrial or/and 
astrophysical origin present causes the residuals, variances, and covariances of certain 
parameters to increase with time. Hence, the low frequency noise corrupts our obser- 
vations and reduces experimental possibilities for better tests of General Relativity 
Theory. At the same time, the rate of growth of residuals and variances of parameters 
can give a valuable information about the red noise itself. We also treat in detail the 
influence of a polynomial drift of noise on the residuals and fitting parameters in order 
to avoid confusion with red noise without the polynomial drift. Results of the analitic 
analysis are used for discussion of a statistic describing stabilities of kinematic (PT) 
and dynamic (BPT) pulsar time scales. 
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2 Sergei M. Kopeikin 

1 INTRODUCTION 

Timing observations of single and, especially, binary millisecond pulsars are widely recognized as being extremely important 
for progressing a number of branches of modern astronomy and time keeping metrology. In particular, implication of pulsar 
timing for testing General Relativity in the strong- field regime (Taylor & Weisberg 1982, 1989), creation of a new astronomical 
time scale based on the high-stable rotation of millisecond pulsars (Ilyasov et al. 1989, Kaspi et al. 1994) as well as setting 
up the upper limit on the energy density of stochastic gravitational wave background in early universe (McHugh et al. 1996, 
Kopeikin 1997a, and references therein) are among the most successful and stimulating recent achivements in this area of 
research. 

The accuracy of pulsar timing observations is now approaching ~ 100 nanoseconds. Such high precision requires con- 
struction of adequate data processing algorithms taking proper account of the relevant physical effects which can contribute 
to the timing signal being processed. Significant theoretical understanding of relativistic celestial mechanics of binary systems 
has been made since the discovery by Hulse and Taylor (1975) of the first binary pulsar PSR B1913+16 (Damour et al. 
1989, Damour & Schafer 1988, Ohta & Kimura 1988, Schafer & Wex 1993, Kopeikin & Potapov 1994, Schafer 1995, Wex 
1995). Presently, one is able to predict both secular evolution and periodic perturbations of the binary system's orbit up to 
the 2^ post-Newtonian approximation (PNA) where emission of gravitational waves starts to influence the orbital dynamics 
(Grishchuk & Kopeikin 1983, Damour 1983). Presently, efforts are being made towards developing theory of motion of binary 
systems in 3 PNA and 3| PNA of General Relativity (Iyer & Will 1995, Jaranowski & Schafer 1997, Jaranowski 1997, Rieth 
1997). 

It is worth remembering that processing of pulsar timing data comprised of topocentric time of arrivals (TOA) of pulsar 
pulses is based on the \ 2 minimization procedure of fitting observational parameters of the pulsar to the adopted model of the 
observed TOA. This includes polar motion corrections UT1-UTC given by the International Earth Rotation Service (IERS), 
the general relativistic model of motion of the observer around the barycenter of the Solar system (Standish 1990) and the 
pulsar around the barycenter of the binary system (Damour et al. 1989), the post-Newtonian transformations between different 
time scales (Brumberg Kopeikin 1990, Fukushima 1995) used in the model, and the law of propagation of electromagnetic 
signals in gravitational fields (Shapiro 1964), as well as the interstellar and interplanetary medium (Rickett 1990, 1996). The 
overall model is rather sophisticated and presently exists in the form of two independently developed computer programs: 
TEMPO (Taylor & Weisberg 1989), and TIMAPR (Doroshenko & Kopeikin 1990, 1995), both being available on the World 
Wide Web. It is worthwhile pointing out that TIMAPR is based on the unique theoretical approach for construction of 
relativistic time scales and reference frames in the Solar system developed in (Kopeikin 1988, Brumberg & Kopeikin 1989, 
1990). 

Usually, the procedure for estimating pulsar parameters is based on the premise that white noise dominates in TOA 
residuals. However, long-term monitoring of certain pulsars has revealed the presence of non-white component of noise of 
astrophysical origin as well (Cordes & Downs 1985, Kaspi et al. 1994, D'Allesandro et al. 1996). It is customary call such 
correlated noise as a coloured Gaussian (or, simply, as "red") one because its spectrum diverges at zero frequency ("the 
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Millisecond and Binary Pulsars as Nature 's Frequency Standards 3 

infra-red catastrophe"). At low frequencies the red noise has a non-flat spectrum and can be described in framework of single- 
or multiple- component power-law model. The lower the timing activity of the pulsar, the further toward low frequencies 
one must look in order to detect the red noise in the spectrum of its TOA residuals. Although being fairly difficult for 
detection, red noise contains invaluable information about diversity of physical processes which take place in the neutron star 
interior, the interstellar medium, early universe, and man-made terrestrial clocks (Lorimer 1996). For this reason, developing 
a rigorous procedure for accounting for the red noise component in TOA residuals is worthwhile. Physicists have been working 
on the problem of adequate treatment of red noise for a long time (see, for example, Stratonovich & Sosulin 1966, Van Trees 
1968, Hooge 1976, Planat et al. 1996). However, one of the first attempts in developing proper statistical procedure of fitting 
pulsar parameters in the presence of low-frequency, red noise has been undertaken only recently by Bertotti et al. (1983), 
and Blandford et al. (1984) in order to obtain unbiased estimates of the pulsar's parameters as well as upper limit on the 
energy density of the GWB radiation. There is a need to further improve on probabilistic models of the red noise and methods 
of estimation of its characteristics containing essential information about long-term processes having both geophysical and 
astrophysical origins, including the stochastic background of primordial gravitational radiation. The present state-of-the art 
of statistical analysis of pulsar timing observational data in the presence of red noise has not yet reached the required level 
of clarity and completeness, and more elaborate methodology has to be developed. This point has been especially stressed by 
Bertotti et at (1983), Hogan & Rees (1984), Blandford et at (1984), and Kopeikin (1997a, 1997b). 

Systematic studies of timing noise in post-fit residuals was started by Groth (1975) and Cordes (1978, 1980) (see also 
Cordes & Helfand (1980), Cordes & Greenstein (1981), and Cordes & Downs (1985)), who put forward a random walk 
noise model along with the use of orthonormal polynomials for analyzing properties of the noise in the time domain. These 
investigations stimulated development of a somewhat equivalent approach for estimating the spectral power of the low- 
frequency noise which was evaluated in the Fourier frequency domain and non-uniformly sampled data by Deeter & Boynton 
(1982) and Deeter (1984). The power spectrum technique was later applied to pulsar timing data by Deeter et al. (1989). This 
formalism was worked further into a practical procedure by Stinebring et al. (1989) and extensively employed by Kaspi et al. 
(1994), and Thorsett & Dewey (1996) for setting an upper limit on the cosmological parameter fi 9 using the Neyman- Pearson 
test of statistical hypothesis. McHugh et al. (1996) recently refined this procedure to make it even more rigorous using the 
approach based on the Bayesian statistic. 

It is worth noting that the analysis in the time domain is more informative than that in the frequency domain. This is 
because noise contains usually both stationary and non-stationary components, and spectral analysis of the noise in frequency 
domain is adequate if and only if the noise itself is stationary (or its increments). For this reason, the procedure advanced 
by Stinebring et al. (1989) is principally restricted by the implicit assumption on stationarity of timing noise and, therefore, 
its implications for processing real observational data are jeopardized. Fortunately, we have been able to prove (Kopeikin 
1997a) and reconfirm in the present paper that under rather general circumstances the procedure of fitting of pulsar spin 
and orbital parameters acts as a low-frequency filter of the non-stationary component of noise which means that pulsar 
post-fit residuals bear, in fact, only the stationary component of noise and the spectral analysis in the Fourier domain is 
© 1999 RAS, MNRAS 000, 



4 Sergei M. Kopeikin 

legitimate. Nevertheless, non-stationary part of noise intrudes into observed values of pulsar spin-down parameters, increases 
their variances, and bringing on large correlations with other parameters. 

Another thing to notice is that inherent to Deeter-Boynton's procedure are certain limitations, which must be clearly 
understood and accounted for in practical implementations. The foremost restriction is that it provides a satisfactory guide to 
the selection as well as construction of practical estimators of a power spectrum of noise as long as the noise is approximately 
represented by a (2r)-th power law; that is, then it is possible to write the noise spectrum down as S(f) = K2 T f~ 21 , where 
r = 1,2,3,..., and any of Kiv is a constant. Let us note that the noise in question can be treated as a random walk in 
rotational phase, frequency, or frequency derivative of the pulsar. Stationary stochastic processes possessing power spectra 
with odd spectral indices S(f) = K.2 T +if~ 2T ~ 1 , (r = 0,1,2,...) are called flicker noise and require development of a more 
advanced statistical approach. A possible way towards the development of such an approach in the frequency domain was 
explored by Blandford et al. (1984). However, it is not as far-going as one would desire. From the mathematical point of 
view the difficulty is due to the existence of an algebraic singularity in the spectrum of low-frequency noise as frequency 
approaches the point / = 0. When the spectral index of the noise is large enough, the singular behavior of the power spectrum 
leads to formally divergent integrals describing timing residuals and variances of measured parameters. To avoid this problem, 
a special regularization procedure must be applied to these integrals to ascribe them definite numerical values and, as a 
consequence, a real physical meaning. There are several known methods of dealing with regularization of functionals with 
algebraic singularities. Blandford et al. (1984) used the simplest procedure of regularization based on the effective spectral 
cutoff of the divergent integrals. As a consequence, their estimators of spectral power of the noise as well as variances of 
pulsar parameters depend on the lower cut-off frequency, /mm, and two numbers, A, and, N* , which are used for minimizing 
root-mean-square error in rotational frequency through the post-fit residuals. The cutoff frequency / m i n and numbers A and 
iV* are not constant and depend on the total span of observations, making temporal behavior of the post-fit residuals and 
variances of parameters rather uncertain. In addition, the simple cutoff of divergent integrals leads to enormous leakage of 
the spectral power from the low-frequency tails of the estimators and as a consequence to the wrong determination of the 
magnitude Km of the noise power spectrum. 

We have suggested tackling the problem of divergent integrals using the theory of generalized functions (Gel'fand & 
Shilov 1964) as the most powerful and simple in dealing with functions having singular spectra. For instance, this theory has 
been successfully applied in Quantum Renormalization Theory (Damour 1975) and for calculations of high-order relativistic 
equations of motion of binary pulsars (Damour 1983a, Schafer 1985). Quite recently the theory of generalized functions 
has been applied by Kopeikin (1997a, 1997b) for development of an adequate treatment of low-frequency timing noises with 
negative integer spectral indices. In particular, we have discovered in these papers for the first time that flicker noise caused by 
the stochastic background of primordial gravitational radiation in early universe leads to the appearance of specific logarithmic 
terms in the autocovariance functions as well as in variances of spin parameters of observed pulsars. This result is re-confirmed 
and considerably extended in the present paper accounting for flicker noise in rotational phase and frequency as well. 

As a final critical remark, let us point out that the problem of proper estimation of variances of orbital parameters of 
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binary pulsars in the presence of low-frequency noise has not been yet discussed in full detail. Therefore, the main goal of 
the present paper is to provide a self-consistent theoretical and numerical study of this problem. Here we use the results of 
our previous work (Kopeikin 1997a) in which the general statistical model of red noise has beeen developed. The plan for 
the present paper is as follows. In the next section we describe an analytical timing model which is used in the subsequent 
discussion. The procedure of estimation of variances of measured parameters in the presence of low-frequency red noise is 
outlined in section 3. A brief description of model of a red noise and its autocovariance function are given in section 4. Section 

5 explains computational aspects of our algorithm. Polynomial drift of timing noise is discussed in section 6. The drift-free 
noise model is employed in section 7 for analytical evaluation of variances of spin and orbital parameters of a binary pulsar 
in the presence of white and the low-frequency noises. Finally, the analytic results are used for studying stability of kinematic 
(PT) and dynamic (BPT) pulsar time scales in section 8. 

To complete the analysis the study of spectral sensitivity of binary pulsars in different frequency bands of the noise 
spectrum has to be done. This requires more work which is currently in progress. These results will be published elsewhere. 

2 TIMING MODEL 

We consider a simplified, but still realistic model of arrival time measurements of pulses from a pulsar in a binary system. It 
is assumed that the orbit is circular, and the pulsar rotates around its own axis with angular frequency v v which slows down 
due to the electromagnetic (or whatever) energy losses. It is also taken into account that the orbital frequency of the binary 
system, rib, and its projected semimajor axis, x, have a secular drift caused by emission of gravitational waves from the binary 
(Peters & Mathews 1963, Peters 1964) bringing about the gravitational radiation reaction force (Damour 1983a, Grishchuk 

6 Kopeikin 1983), radial acceleration (Damour & Taylor 1991, Bell & Bailes 1996), and proper motion of the binary in the 
sky (Kopeikin 1996). 

The moment, T , of emission of the jV-th pulsar's pulse relates to the moment, t, of its arrival, measured at the infinite 
electromagnetic frequency, by the equations (Damour & Taylor 1992, Kopeikin 1994): 



We use the following notations: 

• T - pulsar time scale, 

• t - barycentric time at the barycenter of the Solar system, 

• r* - topocentric time of observer, 

• Ac, A_R0, A^q, A_eq, Asq - clock and astrometric corrections (Taylor & Weisberg 1989, Doroshenko & Kopeikin 1990, 
1995) which one assumes to be known precisely, 



D [T + x sin {n b T + a)] = t + cj> (i) + tj>i (t) 



(1) 



t = t* + A c + A R(S + A nQ + A B0 + A S0 . 



(2) 
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6 Sergei M. Kopeikin 

• D - Doppler factor gradually changing due to the acceleration and proper motion of the binary system in the sky [T] 

• a - initial (constant) orbital phase, 

• rib - orbital frequency (rib = 2-k/Pi,), 

• i - angle of inclination of the orbit to the line of sight, 

• x - projected semimajor axis a p of the pulsar's orbit (a; = a p sini/c), 

• c - speed of light, 

• </>o(i) - the gaussian noise of TOA measuring errors, 

• <f>i(t) - low- frequency gaussian noise caused by the long-term instabilities of terrestrial clocks, effects in propagation of 
radio signals in the interstellar medium, and stochastic background of primordial gravitational waves. 

We also suppose that observations start at the time to = 0. Changing of the starting time of observations from to = to 
to 7^ is achieved by the simple translation of time scale: t i — ^ t — to- It is worth emphasizing that for nearly circular orbits 
the initial orbital phase a = u)q — UbTo where 7o is a fiducial moment of time and u)q is the longitude of periastron. This 
relationship remains valid up to the first order of magnitude in eccentricity and represents an excellent approximation in the 
case of binary pulsars like PSR B1855+09 or PSR J1713+0747. 

The rotational phase of the pulsar is given by the polynomial in time 

N{t) = v p T+±v p T 2 + ±v p T 3 + ± u, T 4 + ± i) p T 5 + v p <j> 2 (T) + 0[T 6 ] , (3) 

where v p , i> v , v p , etc. are pulsar's rotational frequency and its time derivatives all referred to the epoch T — 0, the term 0[T 6 ] 
denotes high order derivatives of the rotational phase, and <j>2 (T) is the intrinsic pulsar timing noise in either rotational phase, 
frequency, or frequency derivative. Solving iteratively equation (^) with respect to T and substituting T for the right hand 
side of equation (^) gives a relationship between two observable quantities Af and t: 

Af(t) = Afo + vt + - v t 2 + - v t :i + i v t 4 (4) 
with 

C (t) = f Mt) + Mt) + Mt), (5) 

where No is the initial rotational phase of the pulsar (No — —vto); v, v, u, ... are the pulsar's rotational frequency and its 
time derivatives at the initial epoch to; the projected semimajor axis of the orbit and its time derivatives at the 

epoch to; f, fib, nb, rib, ... are the pulsar's orbital initial phase, orbital frequency and its time derivatives at the epoch to- All 
non- linear terms of order x 2 ,xe, e 2 , vx, vx, etc. are negligible and have been omitted from m). 



I D = — c =, where Vr and V are correspondingly the relative radial and total velcities of the binary system barycenter with 
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Let us underline that secular variations (time derivatives) of pulsar parameters include not only contributions caused 
by the dissipative physical mechanisms like emission of gravitational waves by the binary pulsar, but also depend on the 
kinematical effects caused by the radial acceleration and proper motion of the pulsar in the sky (Damour & Taylor 1991, Bell 
& Bailes 1996, Kopeikin 1994, 1996). 

Another important point which is usually never stressed, but rather crucial, is the treatment of pulsar timing noise (j>2(T). 
The noise in question contributes to random jumps of the pulsar's phase at the moment T of emission of pulse. However, 
observer at the Earth receives the pulse, and consequently, information about the noise function 4>2 much later after the 
radio pulse has moved across the distance separating the pulsar from observer. Thus, the observer is allowed to think about 
the pulsar timing noise as starting long time ago and the question arises about how to account for this nuisance "memory" 
effect. Actually, the "memory" effect problem is not a specific of the timing noise only. It exists in treatment of observations 
being corrupted by any low frequency noise. As far as we know the only paper, where the "memory" effect in pulsar timing 
data processing has been discussed, is (Kopeikin 1997b). Therein, we have proved that the "memory" effect is negligibly 
small under normal circumstances. Namely for this reason we have deliberately omitted any explicit dependence on time T 
in formulae (0), (H) where the only really important time argument is the Solar system barycentric time t. 



3 PROCEDURE OF ESTIMATION OF PULSAR PARAMETERS 

We assume that all observations of the binary pulsar are of a similar quality and weight. Then one defines the timing residuals 



r 



(t) as a difference between the observed number of the pulse, J\T° s , and the number Af(t,0), predicted on the ground of our 



best guess to the prior unknown parameters of timing model (^), divided by the pulsar's rotational frequency v, that is 

^0) = ^-^ , (6) 

where 9 = {6 a , a = 1, 2, ...k} denotes a set of k measured parameters [k = 14 in the model ([|)] which are shown in Table |l[ It 
is worth noting that hereafter we use for the reason of convinuence the time argument u — n a t. 

If a numerical value of the parameter 6 a coincides with its true physical value 6 a , then the set of residuals would represent 
a physically meaningful noise e{t), i.e. 

r(t,§)=e(t). (7) 

In practice, however, the true values of parameters are not attainable and we deal actually with their least square estimates 
9 a . Therefore, observed residuals are fitted to the expression which is a linear function of corrections to the estimates 9 a of a 
priori unknown true values of parameters 9 a . From a Taylor expansion of the timing model in equation (Q) , and the fact that 
r(t,9) — e(t) one obtains 

14 

r(t, 9") = e(t) - Y^M*> °*) + °0£)> ( 8 ) 

a=l 
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8 Sergei M. Kopeikin 



Table 1. List of the basic functions and parameters used in the fitting procedure. Spin parameters 
8Mo,8v, 8 u,8 v,8 u,8 v , fit rotational motion of the pulsar around its own axis. Keplerian parameters 
8x, 8a, 8ni, fit the Keplerian orbital motion of the pulsar about barycenter of the binary system. Post- 
Kcplcrian parameters 8 x, S x, 8 x ,8 n b ,8 n b fit small observable deviations of the pulsar's orbit from the 
Keplerian motion caused by the effects of General Relativity, radial acceleartion, and proper motion of 
barycenter of the binary system with respect to the observer 



Parameter Fitting Function 



SN 



V-i(t) = 1 



2rl b 

/34 = - i T 4r i/>4(t) = u 3 

^ = i±jNr <fe(i) = « 4 

fa 

(3j — — Sx sin a — Sax cos a ipj{t) — cos u 

/3g — — Sx cos a + Sax sin a i/>g (t) = sin u 

(3g — ^ — S x cos a + Sn b x sin a^ ijjg (£) — u sin u 

010 — (^—Sx sin a — Sn^x cos a^j V'lO (*) = u cos u i 

(3 1 1 — — / _£ x sin (T — 5 Tifo x cos cr J V'llf^) — u ^ cos u 

/3l2 — — -j 1 ( —S x cos o" + 8 rift x sin a \ ' i Pl2 (*) — ^ sinu 

/3i3 — — ^ir ( — S x cos cr + S rift x sin tr) ^13 (*) ~ 14 ^ s i n u 

/3l4 = — ^ [—5 x sin cr — 5 n b x cos cr J */'14 (t) — li 3 cos u 



where the quantities f3 a = S8 a = 9* ~ 6 a are the corrections to the presently unknown true values of parameters, and 
tpa(t,9*) = [§^-~\ _„„ are basic fitting functions of the timing model. 

In the following it is more convenient to regard the increments f3 a as new parameters whose values are to be determined 
from the fitting procedure. The parameters (3 a and fitting functions are summarized in Table [l] with asterisks omitted and 
time t is replaced for convenience by the function u = ntt which is the current value of orbital phase. It is worth emphasizing 
that the basic functions f/>2a-i, (a = 1, 7) are even , and ip2a, a = 1, 7 are odd. We restrict the model to 14 parameters 
since in practice only the first several parameters of the model are significant in fitting to the rotational and orbital phases 
over the available time span of observations. It is also important to understand that the smaller amount of fitting parameters 
one takes, the more significant the contribution of non-stationary part of low frequency noise (see discussion at the end of 
this section). 
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Millisecond and Binary Pulsars as Nature's Frequency Standards 9 

Now suppose that we measure m equally spaced and comparably accurate arrival times each orbit for a total of TV orbital 
revolutions, so we have mN residuals r, = r(ti), i = l,...,mN. Standard least squares procedure (Bard 1974) gives the best 
fitting solution for estimates of the parameters /3 a 

14 mN 

MT) = J2J2 L ^M*Mti), a = l,...,14, (9) 

6=1 i=l 

where the matrix of information is 

mN 

L ab (T) = Y^Mu)Mt*), (10) 

1=1 

the matrix is its inverse, and T — NPt is a total span of observational time. Matrices L ao and are given up to 
numerical factor i n Tables |3[t|]. 

Let the angular brackets denote an ensemble average over many realizations of the observational procedure. Hereafter, 
we assume that the ensemble average of the noise e(t) is equal to zero ^. Hence, the mean value of any parameter j3 a is equal 
to zero as well, i.e. 

< e(t) >= — > < p a >= 0. (11) 

The covariance matrix M aa = < p a Pt > of the parameter estimates is now given by the expression 



M ob (r) = ^^L- 1 L- 1 

c=l d=l 



mN mN 



i=l 3=1 



(12) 



where R(ti,tj) = < e(ti)e(tj) > is the autocovariance function of the stochastic process e(t). The covariance matrix is 
symmetric {M ao = A/j, a ) , elements of its main diagonal give variations (or dispersions) of measured parameters op a = M aa =< 
Pa >i an d the off-diagonal terms represent the degree of statistic covariance (or correlation) between them. Covariance matrices 
are given explicitly in Section 6. 

Subtraction of the adopted model from the observational data leads to the residuals which are dominated by the random 
fluctuations only. An expression for the mean-square residuals after subtracting the best-fitting solution for the estimates Q 
is given by the formula 

mN mN 

< r 2 (T) >= -^J2Y, F ^^ R(U ^' ( 13 ) 

i=i j=i 

where the function 

14 14 

F(ti,t s ) = Sij - 5^£-fctf.(ti)V%(*i), (14) 

a=X 6=1 

is called the filter function (Blandford et al. 1984). These expressions merely demonstrate that the amount of the background 
noise is reduced by the fit for the pulsar's spin and orbital parameters (Bertotti et al. 1983, Blandford et al. 1984, Bard 1974). 

t see Section 6 where the influence of time polynomial drift of the ensemble average of noise is treated in more detail 
© 1999 RAS, MNRAS 000, I 



10 Sergei M. Kopeikin 

Thus, the observed magnitude of residuals is on the average smaller than that of the noise. This is because we have chosen 
the estimates 0* for parameters so as to make the residuals as small as possible. Let us note that the more fitting parameters 
one has in the timing model, the smaller is the mean amplitude of residuals. 

Another remarkable feature of ( |l3| ) is that if the autocovariance function R(ti, tj) contains products of terms of the form 
tpc(ti) x f(tj), where f(tj) is an arbitrary smooth function, they must disappear from the post-fit residuals (Kopeikin 1997a). 
The proof of the statement is based on two exact equalities: 

mN mN mN 

= X>c(*0/(fc), (15) 

i=l j = l i — 1 

and 

mN mN 14 14 mJV 14 14 mN 14 mN 

i=lj=la=lb=l 2=1 a=l b=l j=l 6=1 j=l 

We have set up the hypothesis (Kopeikin 1997a) that the non-stationary part of any low-frequency noise with rational 
power-law spectrum can be represented as a sum of terms being products of a polynomial of time g(t) = ao + ait + a^t 2, + ... 
by a smooth function of time f(t). This has been proved for a fairly general event of low-frequency noise being generated 
by a shot noise random process (Kopeikin 1997b) and seems to be true, if not for all, at least for a considerable number of 
stochastic processes. Taking large enough number of fitting spin-down parameters and corresponding fitting functions ip(t) 
one can represent the non-stationary component of noise as a sum of terms of the form ip(ti) x f(tj). Hence, residuals for the 
example given do not contain the non-stationary component. In real practice, observers fit usually for the first three spin-down 
parameters - initial rotational phase, frequency, and frequency derivative. Such a procedure eliminates from the residuals any 
non- stationary noise components having spectral index n < 6. 

Moreover, we emphasize that the property of the fitting procedure to filter out terms of the form of ipc(ti)f(tj) does not 
actually depend on whether observations are equally spaced. Hence, the conclusion is that, if the autocovariance function has 
a non- stationary component comprised of a sum of terms of the form ip c (ti) x f(tj), one can always make the post-fit residuals 



depending only on the stationary part of the noise 

mN mN 14 14 V mN mN 

i — 1 j — 1 a—1 fa— 1 |_ i= 1 j = l 

For this reason, methods of spectral analysis in frequency domain can be applied without any restriction. 



(17) 



To find the asymptotic behavior of the residuals and elements of the covariance matrix as functions of time (or a number 
of orbital revolutions N) one needs to restrain a model of the stochastic noise process. 

4 MODEL OF NOISE AND ITS AUTOCOVARIANCE FUNCTION 

We have already assumed that noise consists of algebraic sum of mutually uncorrelated components of white and low-frequency 
noises (see Eq.^). White noise is due to measuring uncertainty of TOA while the low- frequency noise has terrestrial and/or 
astrophysical origin and becomes significant only on long enough span of observational time. Investigation of its nature is 
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one of the most important problems of pulsar timing since it gives us a key to much better understanding of the very deep 
fundamentals of physical laws governing the evolution of Nature on long time intervals. Especially important in this respect 
are observations of millisecond pulsars in binary systems having the best attainable accuracy and, consequently, the low level 
of white noise component in TOA residuals. 

White noise is stationary and has a flat constant spectrum 

S(f) = ho, (18) 

with ho being a fixed parameter determined from the measured level of white noise on short time scales (see Eq. (^). The 
autocovariance function of white noise is 

R(T) = h S(r), (19) 

where 8(t) is the Dirac delta function and r = ti — tj with U, tj being moments of TOA of the i-th and j-th observations. 

The generalized model of low frequency noise was developed by Kopeikin (1997b) using the shot-noise approximation. 
It can describe both random walk and flicker noises having, as a rule, different physical origins. The autocovariance function 
of the noise includes both stationary and non-stationary components which are treated on equal footing. We suppose that 
the autocovariance function of any stochastic process under consideration is a function of real variable and can be split 
algebraically into stationary and non-stationary components 

R(U,tj) = R + (U,tj)+R-(ti,tj). (20) 

The function R + (ti,tj) describes the non-stationary part of the noise and R~(tt,tj) = R~(\ti — tj\) is its stationary counter- 
part. The latter can be displayed in the frequency domain through the cosine Fourier transform 



R~{r) = 2 S(f) cos(27r/r)df, T = U-tj (21) 

Jo 

where S(f) is the spectrum of R~{t) 

f oo 

S(f) = 2 R-(r) cos(2ty fr)dT. (22) 
Jo 

The stationary part of the red noise has a power-law spectrum, S(f), being proportional to f~ n where n is called the spectral 
index of the noise and in the present paper n = 1, 2, 3, 4, 5, 6. Corresponding autocovariance functions of the stationary com- 
ponents are represented by either polynomial of time (for random walk noise) or polynomial of time plus logarithmic function 
of time multiplied by another polynomial of time (for flicker noise). Nomenclature of the corresponding names of noises, 
their spectra, and stationary part of the autocovariance functions are given in Table ^ which is an extract from (Kopeikin 
1997b). Non-stationary parts of the autocovariance functions are expressed as a sum of terms of the form: [basic function 
of fitting procedure ipc(ti[) x [a smooth function f{tj)\ and are given below in Section 7. As mentioned above, such terms 
do not contribute to post-fit timing residuals. However, they do contribute to the fitted values of the pulsar's spin parameters. 
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Table 2 . The spectra of timing noise and the stationary part of the autocovariance functions R~(t{, tj ) 
used in the paper. Constant parameters h n , where n = 1,2, ...,6, characterize the magnitude of the noise 
spectrum. The quantity r = tj — tj , and one introduces notation r_ = — r, if r < 0, and t_ = 0, if r > 0. 



Noise Spectrum Autocovariance Function 



White Noise (WN) 


h 




h S(T) 


Flicker Noise in Phase (FPN) 


hlf 


-l 




Random Walk in Phase (RWPN) 


h 2 f 


-2 


— h 2 T_ 


Flicker Noise in Frequency (FFN) 


h 3 f 


-3 


£-h 3 T 2 ln\r\ 


Random Walk in Frequency (RWFN) 


hif 


-4 




Flicker Noise in Frequency Derivative (FSN) 


h 5 f 


-5 


-sfcr^ 7 " 4 ln l r l 


Random Walk in Frequency Derivative (RWSX ) 


h s f 


-6 


-T5o h <3 T - 



Before going further, let us discuss in more detail computational aspects of our analytical approach for calculation of 
residuals and covariance matrices under the influence of different noise processes. 

5 COMPUTATIONAL ASPECTS 

In order to accomplish all calculations analytically we assume that the number of observations is so large that any sum over 
observing points can be approximated by the integral over the observing period T* . For calculational convenience, we assume 
that observations are commenced at the moment to = — T*/2 and finished at to = T*/2 so that the interval of integration 
time is symmetric with respect to the origin of time scale. Such a shift of the time scale is always possible and preserves the 
invariance of timing formula (^). Moreover, we assume that observations are equally spaced with small enough interval of time 
between successive observations At = Pb/m where m is the number of observations per one orbital period. A total number 
of orbital revolutions N is supposed to be large (N > 30). 

Under the given conditions one may apply to any smooth function /(t) 6 C°°[— to, to] the Euler-Maclaurin summation 
formula (2.9.15) given in textbook of Davis & Rabinowitz (1984) and having the form of an asymptotic expansion 

mN mN 

^/(^)=^/Ho + (i-l)At] = |l / f(u)du + E(m,N), (23) 

where T = irN, and the equalities ntAt/Pt = ribdu/ '(27r) and T = ri(,to = ntT* /2 have been used^. Herein E(m,N) is the 
error of approximation of the summation formula 



E(m,N) = \ [/(-T) - /(T)] (If)"" 1 [/^(T) - f , (24) 

where derivatives /'■ fe ' of the function /(u) are taken with respect to the argument u, and the constants are called the 
Bernoulli numbers with the numerical values 

§ New number T counts the span of observational time in terms of orbital revolutions of the binary pulsar in question. When discussing 
observations of a single pulsar we treat the number N as being counted in years 
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B 2 = l 



Ba — — — , 



Be, — — , 



Bg = ——, 



etc. 



(25) 



6 30 42' 30 

It is not difficult to check that the error E(m,N) is decreasing, at least as (mN)' 1 comparativly to the main term, or even 
faster. It may be made small enough (less than 5observations per orbital period m > 10 and the number of orbital revolutions 
N > 30) (Wex 1997). We assume, in what follows, that this condition is fulfilled and the residual terms E(m, N) are negligibly 
small. 

Because the autocovariance function consists of a linear combination of two terms representing non-stationary and sta- 
tionary components, the covariance matrix (jl2|) is split into two liearly independent terms as well: 



M ab = M++M ab , 



(26) 



where Mi depends only on R + (ti,tj) and M~ b depends only on R~(ti,tj). Let us underline that the shift in the origin of 
initial epoch of observations leads to the shift of arguments in the autocovariance functions ti ^ ti + T*/2, tj — » tj + T* /2. It 
leaves the stationary part of the autocovariance function R~(ti,tj) invariant, depending on the difference ti — tj only, though 
produces changes in the non-stationary part of R + (ti,tj) (see section 7 for more detail). 

Calculations of integrals depending on the non-stationary part R + (ti,tj) of the autocovariance function are relatively easy 
to handle. Let us remind that the function R + (ti,tj) for any noise model is composed of the finite sum of products of basic 
functions ifii(ti)ipk(tj), (J, k < 6) in the event of random walk noise to which functions having the structure tpi(ti) \n(n b tj +T), 
/ < 3 are added in the event of flicker noise (Kopeikin 1997b). On one hand, for random walk the corresponding part of the 
covariance matrix will be 



(27) 



SalSb 



(J,ft<6), 



where 5 ab is the unit matrix. On the other hand, in the event of flicker noise one has additional contributions 



M ab ~ 



J2Mh)Hn b t 3 +T) 



3=1 



(28) 



Sal 



14 T 

%J2 L ™ Vto(u)ln(« + T)du 

A— 1 J —T 



(I < 3). 

Relationships (p^)-^3) lead to the important conclusion that, for the timing model under consideration, the non-stationary 
part of any autocovariance function contributes only to the elements of the covariance matrix M a b with indices a, b < 6 for the 
random walk noise, and to the elements of M a b = Mba with a < 3, b = 1, 2, ...14 for the flicker noise. This specific behaviour 
of non- stationary part of autocovariance functions was noted partially by Groth (1975) in the context of his investigation of 
random walk processes in pulsar timing based upon the usage of orthogonal polynomials. Our derivation of the given result is, 
in fact, more general as it does not depend on the particular choice of special functions used for expansion of random process 
and is applicable both to random walk and flicker noises as well. 

Calculations of integrals depending only on the stationary part of autocovariance function are fulfilled after making a 
© 1999 RAS, MNRAS 000, I 
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Figure 1. Geometrical illustration showing change of independent variables and transformation of the domain of integration in the phase 
space of two time arguments. Dashed arrowed lines indicate paths of integration. 

preliminary transformation of independent variables (see Fig. Q). Let us denote u = ribti, v = n^tj, t% 7^ tj and make a 
transformation 

u + v 



V = 



(29) 



2 ' 2 

It is obvious that the stationary part R~(ti,tj) of the autocovariance function depends only on \x\. Hence, the part of the 



covariance matrix M ab is transformed after using formula (E9I) to the integral 



4> a {u)tpb(v)R (u — v)dudv = / dxR (x) / dyA a b(x,y), 
■i Jo Jo 



where the symmetric matrix 

A ab (x, y) = -2 [1 + (-l) a+b ] [ip a (u)i> b {v) + 



(30) 



(31) 



The result shows that in the analytical approximation under consideration the elements of matrix M~ b 7^ 0, if and only if the 
basic functions %p a {u) and t[>b(v) are both either odd or even. In other words, any fitting parameter having odd number has 
no correlation with that having even number. 

For this reason, when calculating integrals by formula ( ^l| ) it is more convenient to re-group fitting parameters in order 
to represent the matrix A a b{x,y) into a block-diagonal form 



A ab 



A! 
A 2 



(32) 



where the left upper block A\ — A2a-l,2b-i depends only on odd fitting functions, and the right lower block A2 = Ai a ^.b 
depends only on even ones (a, b — 1, 7. Due to this re-arrangement of parameters, the matrix L a b and consequently the 
inverse matrix L~ b are also reduced to the block-diagonal form having the same structure as in (^). The advantage of this 
re-arrangement is that one can work with each of the block matrices independently, simplifying the calculations considerably. 

© 1999 RAS, MNRAS 000, §4^ 



Millisecond and Binary Pulsars as Nature's Frequency Standards 15 

Nevertheless, it is worthwhile to remind that such simplification can be achieved in the main approximation only. If one took 
terms of higher order in the expansion ( p3| ) the re-arrangement of elements of matrix M~ b would be worthless since odd and 
even fitting parameters get correlated. In such event one can make a progress in calculations only if numerical methods are 
implied. This is the case one meets in a real practice. 

In the analytical approximation under consideration the integral (|3(]) was calulated using the enhanced version 2.2.3 of 
MATHEMATICA for Microsoft Windows. It turns out that the outcome of calculation can be always expressed either as a 
polynomial and trigonometric functions of time in the event of random walk, or as a polynomial, trigenometric functions, and 
sine/cosine-integrals in the event of flicker noise. Since the number of observations mN has been assumed to be very large 
all functions appearing in elements of covariance matrix ( |3o| have been expanded into asymptotic series with respect to the 
small parameter e = ^ = ^ . Then, only the first term of the expansions has been retained and all residual terms have been 
abandoned. An explicit result of such calcualtion is, for instance, the matrix of information L a b = j^C a b given in Tables ^ 
and |^. Analytical expressions for residuals and elements of covariance matrices of pulsar's fitting parameters in presence of 
white and red noises are given in Section 7. 



6 TREATMENT OF POLYNOMIAL DRIFT OF TIMING NOISE 

For mathematical convenience, in this section we shall assume that the number of spin-down parameters of timing model 

is equal to M and the total number of fitting parameters is A" p]. Up to now, we have been concerned mainly with the 

mathematical formulation of procedure of estimation of pulsar's parameters under the assumption that noise has only random 

fluctuations, that is equation ([n]) is true. However, nonrandom variations (drifts) of the noise do exist and can be sometimes 

modeled by deterministic polynomial functions of time as it takes place in description of electromagnetic breaking of pulsar's 

rotation. Hence, it is conceivable to imagine that the ensemble average of noise is modeled by a (P — l)-th degree polynomial 
p M 

< e(t) >= '^^ct k t k ' 1 = ^2, n b k+1 a k^k{t) + {o-M+itf + ... + aptf _1 ) , (33) 

k=l k=l 

where oik are constants which numerical values get smaller as k increases. From the definition of fitting parameters and fitting 
functions as well as equation ^ it is easy to verify that the mean value of the parameters reads as 

K P mN K M mN K P mN 

< p a {T) > = Y£22 L * ak M t *rt~ 1 = EEE^ 1 ^ 1 ^^^^ +EE 5>«6W*(tot? -1 

b=lk=li=l b=lk=li=l b=l h=M+l i=l 

(34) 

M K mN 

= J2 n b k+ls *ka k + Y^Z L -bMU) (a M+ itf + ... + aptf" 1 ) . 

k=l 6=1 i=l 

Thus, we conclude that until number P of the polynomial's coefficients (^) is less or equal than number M of spin-down 
parameters being used in deterministic timing model it changes only the mean values of the first M parameters. In such event 
one can make several more iterations in the procedure of minimization of square of residuals in order to eliminate completely 



' In particular, the numbers M and K are equal to M 
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the polynomial drift from the ensemble average of the noise by means of re-shifting first P spin-down parameters. As a result, 
one will get < r(t) >= 0. 
If P > M, then 



< r{t) >= a M +it M + ... + a P t p 1 - ^^L~ 6 Va(t) 



E>(**)( 



,M 

CXM+ltj + . 



(35) 



One recognizes that in this event < r(t) >^ though in practice it can be already within the limit of rounding error of 
numerical computations. It is worth emphasizing that even if it is the case, some fitted parameters f3 a can, nevertheless, have 
the mean value lying above the limit of the rounding error. We conclude that it is impossible to make the mean value for the 
residuals and all fitting parameters equals to zero if the number of spin-down parameters in timing model is less than number 
of coefficients in polynomial drift of the ensemble average of noise. This conclusion holds true for both white and red noise 
regimes. 

Let us turn attention now to the calculation of residuals in the presence of polynomial drift of the noise's ensemble average 

(p3[). The autocovariance function of noise reads now as 

p p 

< e(t ( )e(tj) >= R{U,tj) +J2J2 a k<xitt 1 t l f\ (36) 

k=l 1=1 

where R(ti,tj) is the autocovariance function of the noise with zero- valued ensemble average. Influnce of the function R(ti,tj) 
on residuals and covariance matrix of fitting parameters is explicitly demonstrated in the next section. Calculation of ensemble 
average of residuals with making use of expression ( |36| ) and the property that fitting procedure filters out all terms of the 
form i> a {U)f{tj) yields 

mN mN 



<^)>=^EE^. 



i=l 2=1 



R {ti,tj 



E E akait ^ ^ 1 



k=M+ll=M+l 



(37) 



Thus, one can see that if the condition P < M + 1 holds, then the polynomial drift of noise does not contribute to the 
ensemble mean value of square of timing residuals. In the opposite situation, when P > M + 1, there extra contribution exist 
and there is a danger of confusion of white noise having a polynomial drift with a red noise without such a drift. 

Regarding covariance matrix of parameters in presence of polynomial drfit of noise it is not so difficult to generalization 
( p^ ) using expression (^6|) for the autocovariance function of noise. It results in 



M ab = M ab + Ml + AM ab 



(38) 



where the stationary part, M ab , the non- stationary part, M^ b , and the increment AM a b of covariance matrix read, corre- 



spondingly, as 

K K 

M-(T)=Y,Y, L ^ L ^ 



c=l d=l 

K K 



»=i j=i 

mN mN M M 



(39) 



i=l j = l 



+2 



-Hit „ 

, OpkXXk 



EE L ^^) E ait 



i-i 

3 



c=l 3 = 1 



l=M+l 



(40) 
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(41) 



'Y^jii>c(t i )ii) d (tj) ^2 ^2 afc«i** 1 t l j 

c=ld=l Li=lj = l fc=M+i;=A/+l 

Herein, we have assumed that the non-stationary part of autocovariance function of noise R + (ti,tj) consists only of products 
i>a{ti)f{tj) with a < M. 

Analayzing the structure of relationships (|3^)-(|4l|) we draw the following conclusions: 1) if the number of coefficients P 
in the polynomial drift of noise is less than or equal to the number M of spin-down parameters in the timing model, then 
it changes only non-stationary part M* b of the covariance matrix for spin-down parameters; -2) if the number of coefficients 
P in the polynomial drift of noise is more than or equal to the number M of spin-down parameters in the timing model, 
then it influences the non-stationary part of the covariance matrix for spin-down and orbital parameters and gives rise to 
AM„t. Hence, if we do not account for the polynomial drift properly it may worsen numerical boundaries for true estimates 
of variances and covariances of fitting parameters. 

At this point, it is worthwhile to emphasize that in order to make effects of polynomial drift of noise as small as possible 
one has to introduce to the timing model as many spin-down parameters as it is required by observational accuracy and 
optimization of the least-square minimization procedure. 



7 RESIDUALS, COVARIANCE MATRIX, AND VARIANCES OF MEASURED PARAMETERS 
7.1 General Comments 

In what follows we assume that effect of fitting spin-down parameters is so suppressive that the polynomial drift of the noise is 
completely eliminated from all residuals and mean values of parameters. When calculating covariance matrices of parameters, 
it is worth keeping in mind that if one takes only stationary components of low frequency noise it is not sufficient to demand the 
variances of all fitting parameters to be positive since the first several spin-down parameters can have variances with negative 
values (see Tables (9)-(16) for more detail). Negative variances are, of course, physically meaningless. Hence, we conclude that 
it is incorrect to disregard the non-stationary part of red noise for it hampers the physical interpretation of these variances. It 
turns out that, in order to obtain positive values for variances of all fitted parameters, one must account for the contribution 
from the non-stationary component of noise as well. This clearly demonstrates the role of the non-stationary noise component 
in the fitting procedure. The energy distribution amongst the fitted parameters is such that the non- stationary part of red 
noise contains as large amount of energy as the stationary part does (or even more) for the first several parameters. 

An additional point to note is that, if the number of spin-down parameters is not big enough, the non-stationary part 
of the autocovariance function can not be represented as a linear combination of products of a smooth function f(t) by the 
fitting function t[> a (t), and complete filtering out non-stationary component of noise from residuals is impossible. In such a 
case, if only the stationary part of the noise is accounted for and the non-stationary part is omitted, a negative mean value 
of residuals can be obtained that is physically unadmissible. Taking into account the non- stationary noise component brings 
the mean value of residuals back to the positive value. By this consideration, we would like to underline the role which 
non- stationary part of noise plays in fitting procedure. We consider those models of low frequency noise which take into 
© 1999 RAS, MNRAS 000, 
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Table 3. Elements of matrix of information C a f,. Quantities T = ttN and Q = cosT. 





1 


3 


5 


7 


9 


11 13 


1 


2T 












3 


2T 3 
3 


2T 5 










5 


2T 5 


2T 7 


2T 9 








7 





4TQ 


8T 3 Q 


T 






9 


-2TQ 


-2T 3 Q 


-2T 5 Q 


T 

2 


^2~ 




11 


4TQ 


8T 3 Q 


12T 5 Q 




2~ 




13 


-2T 3 Q 
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-2T 7 Q 


2~ 




2~ ~7~ 



account only stationary component to be incomplete and can lead to erroneous results. For this reason, we are very careful 
in dealing with the non-stationary noise component in order to properly take into account or to prove that it is insignificant 
for calculations. In this context, let us note once again that we have made a shift of the time scale which changes the orbital 
phase u i— > u + ttT and introduces additional terms in the non-stationary part of the autocovariance function in comparison 
with expressions for R + (ti,tj) having been given in Kopeikin (1997b). 

With these comments we are ready to demonstrate how white and red noise affect pulsar timing residuals and covariance 
matrices of fitting parameters. 

7.2 White Phase Noise 

The spectral power of white noise is constant. For this reason, the timing residuals are constant as well 

< r 2 >= h , (42) 

and the covariance matrix of measured parameters M a ;, coincides with the inverse matrix of information ( Jic| ) L~£ . The matrix 
of information L a t — ^C a b, and the elements of the inverse matrix are given in Tables ^ - ^[ 

On multiplying the main diagonal of by the magnitude of the measurement error ho yields the variances of the measured 
parameters. One can see that precision of the parameter determination depends both on the rate of making observations m and 
on the number of orbital revolutions iV during the observational session. This is a distinctive feature of timing measurements 
in the event when white noise is the only source of errors. We also note that as the observational span T increses the variances 
of all parameters decrease. Thus, until white noise dominates in observed timing residuals, the accuracy of determination of all 
parameters will be improved. It is worth emphasizing as well that there is a similarity in the precision of determination of spin 
and orbital parameters. Indeed, if one compares variances for the rotational phase and its time derivatives with corresponding 
quantities for orbital motion one finds that dependence on T is essentially the same. For making this statement more obvious 
we give relationships between orbital parameters and parameters of the timing model shown in Table 1 

< (8a) 2 > = — (< $ > cos 2 a+ < 01 > sin 2 a) , (43) 
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Table 4. Elements of matrix of information C a t- Quantities T = ttN and Q = cosT. 
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Table 5. Elements of the covariance matrix C ^ of pulsar's parameters for white noise in phase. Quantities 
T = ttN and Q = cosT. 
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Table 6. Elements of the 
T = ttN and Q = cosT. 
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<{5n i > > = \ (< Pl > cos 2 a+ < Pi > sin 2 a) , (44) 
n b x 

< {5n f > = ^ (< A 2 i > cos 2 a+ < A 2 2 > sin 2 a) , (45) 



< (5ri b f > 36n£ , 2 ^ 2 . . „2 . . 2 \ 



< (5s) 2 > = (</?!> cos 2 a+ < /3 2 > sin 2 cr) , (47) 

< (5s) 2 > = n 2 (< Pi > cos 2 cr+ < Pfo > sin 2 cr) , (48) 

< (6'x) 2 > = 4ra£ (< /3 2 2 > cos 2 cr+ < ^ > sin 2 a) , (49) 



< (<5'af) 2 > = 36ng (< /3 2 3 > cos 2 a+ < p' 2 4 > sin 2 cr) . (50) 



Calculation of variances of orbital parameters is accomplished using results shown in Tables |3j-|6| It is remarkable that there 
is a symmetry between variances of odd and even parameters of the timing model, so that in case of white noise dominance 
one has, for example, < p 2 > = < /3| >, and so on. Because of this symmetry there is no dependence of variances of orbital 
parameters on the initial orbital phase a. From equations (|43|)-(|5C|) we have 



<(^) 2 > - ^ (51) 



9h 



< {5n b ) 2 > _ 75/io 



< (Sh b ) 2 > 45/ipw 2 
n 2 " x 2 T 5 



< (Sh b f > 1575fe n^ 
n 2 ~ x 2 T 7 



(52) 



(53) 



(54) 



<(^) 2 > = (55) 



9ho 



2 



< [Sx) > = 4T 3 , (56) 



© 1999 RAS, MNRAS 000, 1-^ 



< (Sx) 2 > = 



45/ion* 
T 5 



Millisecond and Binary Pulsars as Nature's Frequency Standards 21 

(57) 



= (58) 

which can be easily compared with variances of corresponding spin-down parameters in Tables ^, ^. In particular, we find 
that ratio of root squares of variances of the spin-down parameters to corresponding parameters describing time evolution 
of orbital phase is approximately equal to x/Pb, that is < 8N$/v > / < STq >~ x/Pt, < 8vjv > / < 5Pt/Pt >~ x/Pb, 
< 8i> /v > / < SPt/Pb >~ x/Pb, etc. As a concequence of this observation we conclude that in the case of domination of white 
noise spin-down parameters of binary pulsars are always determined better than corresponding parameters of orbital phase. 

7.3 Flicker Noise in Phase (1// noise) 

Pulsar timing residuals grow with time according to the relationship 

(59) 



< r 2 >= hi 



where hi is a constant number characterizing intensity of the noise. 

The non-stationary part of autocovariance function for PN random walk is represented as (Kopeikin 1997b): 

R + {U, tj) = hi [ln(u + T) + hi(v + T)] , (60) 

and can be expressed through the basic functions 

R+iU,^) =fti[^i(v)ln(u + T)+V2(u)ln(u + T)]. (61) 

From equation ( ps| ) , the contribution from the non-stationary part of the flicker noise to the covariance matrix is given by the 
integral 

14 „ T 14 »T 

<^2 L bc ip c (u)Mu + T)du + 5ibJ2L-i ip c (u)la(u + T)du . (62) 
It is clear from (^) that the only non-zero elements of the matrix M^ b can be M^~ b = M b \ (b = 1, 14). The most interesting 



Ml = hi 



component is the contribution of the non-stationary part to the variance of the first measured parameter 

M+=-gi+21n(2T) (63) 

The matrix M~ b is generated by the stationary part of the autocovariance function and is given in Tables [B| and Let us note 
that, as expected, variances of all measured parameters including the first one are positive quantities. 

Comparing results of the variance calculations with those obtained in previous section for white noise we note that flicker 
noise in phase worsen our ability in determining variances of spin-down parameters. At the same time flicker noise does 
not disturb variances of orbital parameters and symmetry between pairs of odd and even parameters. For this reason, using 
general formulae (^3[)-(|5C|) one can see that variances of the orbital parameters do not depend on the initial orbital phase a 
and have the same dependence on the total span of observation T as it was in case of white noise. We would like to emphasize 
the appearance of logarithmic terms in timing residuals and the covariance matrices given in Tables ^ and ^. Logarithmic 
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Table 7. Elements of the co-variance matrix M ah of pulsar's parameters for flicker noise in phase. 
Quantities T = nN , Q = cos T, and T = 7 + ln(2T). The magnitude of the noise hi is omitted. 
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4T 3 
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- 105 Q [17 + 6or] 
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6 1T X 
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16T y 



4T 5 



. 3157T 

4T 5 



10571 
4T 5 



175j 
4T 7 



Hereafter the letter 7 = 0.577215... denotes the Euler's constant. 



Table 8. Elements of the covariance matrix M ab of pulsar's parameters for flicker noise in phase. Quantities T = 7riV, 
Q = cos T, and T = 7 + ln(2T). The magnitude of noise hi is omitted. 



_ 3381 
64T 4 

46893 
1280T 

_ 693Q 
160T2 

315Q 



[61 + 20r] 



nrylQ 
32T 4 



r [261 + ioor] 



I 2 

64T 6 

- 4851 
32T« 

357Q 
16T 4 
_ 105Q 



[341 4- 180r] 



1785Q 



7:i5U 
: : ;2T S 



[91 + 60r] 



161)083 
1280T 1U 

_ 693Q 
32T b 



2079O 



[21 + 2or] 



3465Q 



■i|iS[43 + 60r] 
64Tl () 



r [289 - 120r] 



-tit 

2T5 



. in5 7r 

4T 5 



- [49 — 20r] 



terms are characteristic for any flicker noise as it will be shown in subsequent sections. However, practical observation of the 
logarithmic terms in temporal behaviour of residuals and variances of measured parameters is a challenge for observers since 
it is difficult to distinguish InT from constant if observational span is not long enough. Perhaps, this explains why logarithmic 
behavior of timing residuals and variances have not been yet found and flicker noise with spectral index s = 1 has not been 
ever identified in pulsar timing observations. 

7.4 Random Walk in Phase (l// 2 noise) 



The timing residuals grow proportionally to the number of orbital revolutions iV 
< >= ^^T, 



(64) 



where /12 is a constant number characterizing intensity of the noise. The non-stationary part of autocovariance function for 
PN random walk is represented as (Kopeikin 1997b): 
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Table 9. Elements of the covariance matrix M , of pulsar's parameters for random walk in phase. Quantities 
T = ttN and Q = cos T. The magnitude of noise h 2 is omitted. 
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Table 10. Elements of the covariance matrix M ab of pulsar's parameters for random walk in phase. Quan- 
tities T = ttN and Q = cos T. The magnitude of noise /12 is omitted. 
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R+iu,^) = h 2 T {vi(^ife) + ^ [M^)Mt]) + Mu)Mt])]} ■ 

Using Eq. ifoW) the covariance matrix for non-stationary part of PN results in 



1 1 

SalSbl + ■^5al8b2 + 7^SblSa2 



(65) 



(66) 



Stationary part of the PN noise gives the covariance matrix displayed in Tables [)] and [l^. One can see that the random 
walk in phase makes it impossible to track the rotational phase of the pulsar after amount of orbital revolutions exeeds 
(28 167r)/(659// 2 h 2 ) ~ 13A/{v 2 h 2 ) since the moment when random walk in phase becomes the dominant source of noise in 
timing residuals. In this case the parameter A/"o becomes non-informative|JJj All other spin-down parameters can be measured 
but accuracy of their measurement is lower comparatively to the cases of white and/or flicker noises in phase. On the other 



II A parameter is called informative if its mean numerical value obtained in the fitting procedure is much less than its variance. In the 
opposite case, the parameter is called non-informative. We also call the initial rotational and orbital phases, Wo and <r, informative 
parameters if their variances do not exeed 2tt. From a physical point of view, some of the fitting parameters become non-informative 
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hand, time dependence of variances of orbital parameters is the same as it was in the case of white noise. Hence, we are 
still able to determine the initial orbital phase a, orbital frequency rib, and so on. However, the symmetry between pairs of 
corresponding odd and even parameters of the timing model under consideration is not preserved. For this reason, accuracy 
in determination of variances of orbital parameters slightly depends on the initial orbital phase (which remains to be the 
informative parameter) . 



7.5 Flicker Noise in Frequency (l// 3 noise) 



The timing residuals are proportional to the square of the number of orbital revolutions 
<r2> =W5 h3T2 > 



(67) 



where /13 is a constant number characterizing intensity of the noise. The non-stationary part of autocovariance function for 
flicker noise in frequency is represented as (Kopeikin 1997b): 



R + {U,tj) = |/i 3 T 2 |-^i(«)^i(w) - ±i/}i(v)ih(u) + -±ih(v)Mu>) ~ 4?Mv)i>2(u)+ 

[tpi(v) + |Va(«) + %ipi(v)u - ^ipi(v)u 2 ] ln(w + T)+ 

+ [Vi(u) + ^i>2{u) + (u)v ~ ^ipi(u)v 2 ] ln(u + T) 
Using Eq. (feq), the covariance matrix for the non-stationary part of the flicker noise in phase results in 



(68) 



= feT 2 ^ —^SalSbl — ^f8 a l5b2 — ■£f8bl8a2 ~ <5 a 2<?>i,2 + 



+ - 



Sai^L^ 1 ip c (u)(l- ^yn(u + T)du + S bl J2Lac Mu) (l- ^)ln(u + T)o 

c=l c=l 

14 pT 14 pT ~| 

Sa2^2L- c 1 Vc(w)(l+^)ln(u + T)dM + 5 i)2 ^L- 1 / ip c (u) (l + ^) In (w + T) du I 

o=l c=l J - T ' 



(69) 



It is clear from that the only non-zero elements of the matrix M^ b are M^ b = M^, M^ b = M^, (b = 1,..., 14). In 
particular, the elements and M^t, which contribute to variances first and second parameters are given 
1991" 



M+ = h 3 T 2 



M.7. 



In (2T) - 
2 In (2T) + 



1680 
143 



1680 



(70) 
(71) 



The stationary part of the covariance matrix M~ b is given in Tables ( |l"l| ) an d ( |l^ ) . Sum of matrices M+ b and M~ b gives positive 
numerical values to variances of all measured parameters. 

Flicker noise in frequency, if it dominates in timing residuals, will not improve the accuracy in determination of the initial 
rotational phase and frequency. Moreover, variance of the initial rotational phase grows as span of observations increases. 
Information about numerical value of the initial rotational phase will be completely lost when the root square of the variance 

when the ratio of red timing noise to the deterministic signal (Q) cxeeds a certain level. The stage at which this happens may be found 
from comparision of the mean values of the parameters and their variances. 
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Table 11. Elements of the covariance matrix M , of pulsar's parameters for flicker noise in frequency. 
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Table 12. Elements of the covariance matrix M . of pulsar's parameters for flicker noise in frequency. 
Quantities T = nN and Q = cos T. The magnitude of noise /13 is omitted. 
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becomes equal to 27r. Taking into account the non-stationary component (|7^) of the covariance matrix we find that it happens 
after amount of orbital revolutions iV exeeds ~ 2/ (vh}J 2 ). For single pulsars, this corresponds to a time interval of about 
~ 2/(uh^ 2 ) years after which the initial rotational phase becomes the non-informative parameter. We also note that variances 
of the rotational frequency v and initial orbital phase a can not be determined better than their values obtained before the 
flicker noise in frequency has commenced to corrupt observations. Symmetry in variances of corresponding pairs of odd and 
even orbital parameters is completely destroyed. One interesting consequence of this phenomena is that the measurement 
precision of orbital parameters in the presence of flicker noise crucially depends on numerical value of the initial orbital 
phase until it remains to be informative parameter. For example, from equations ( ^ ) and ( p2[ ) one can see that if we arrange 
observational data to have a ~ then parameter x is determined better than the orbital frequency rib. Inversly, if observational 
data are prepared to make a ~ 7r/2 then the rotational frequency may be measured better than x. The situation looks similar 
to that existing in the quantum mechanical non-demolition experiments (Braginsky & Khalili 1992) in which the quantum 
© 1999 RAS, MNRAS 000, I 



26 Sergei M. Kopeikin 



Table 13. Elements of the covariance matrix M , of pulsar's parameters for random walk in frequency. 
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system is prepared in such a way to make the quantum measurement of one of the conjugate variables as good as possible at 
the same time losing information about the other conjugate one. 

Certain cases in D'Alessandro et al. (1997), for which the derived spectra of the observed timing residuals show spectral 
indices n = 3,5, could well be indicators of the flicker noise of the respective kind. 

7.6 Random Walk in Frequency (l// 4 noise) 

The timing residuals are proportional to the cube of the number of orbital revolutions 
2 



< r >= 



6435 



h 4 T 3 



(72) 



where is the constant quantity characterizing intensity of the noise. The non-stationary part of the autocovariance function 
for PN random walk is represented (Kopeikin 1997b) as 



(73) 



+ ^1p2(u)lp2(v) - jij [lpl(u)lj)4,(v) + 1p4,(u)lpl(v)] + 1^3 (u)lp2 («) + tp2 (u)lp 3 («)] 

Using (^) the covariance matrix for non-stationary part of PN results in 



Ml = -h 4 T* 



3 3 3 1 1 3 3 

SalSbl + — S al S b2 + — S bl S a2 + 7 ^5 a2 5b2 - -p^S a l5b4 - -p^SblS a 4 - -p^5 a 25b3 + -77p3^2<5 Q 3 



(74) 



The stationary part of the noise gives the covariance matrix displayed in Tables |l^ and jl4| 



Random walk in frequency significantly contributes to variances of the first three spin-down parameters so that the 
initial rotational phase and frequency become non-informative parameters. Rotational frequency derivative continues to be 
the informative parameter but precision of its measurement remains constant and can not be improved anymore. This has a 
direct consequence on our ability to use rotational motion of the pulsar as a time scale (for more detail see Section 8). 
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Table 14. Elements of the covariance matrix M~ h of pulsar's parameters for random walk in frequency. 
Quantities T = ttN and Q = cos T. The magnitude of noise /14 is omitted. 
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Random walk in frequency also converts the initial orbital phase or semi-major axis of the binary pulsar from the 
informative to the non-informative parameters. Indeed, equations (|43|), ( ^?| ) and variances of parameters P7,f3s taken from 
Tables yield 



< (da)' > 



hi / 2493 ___ 2 _ , 3T : . 



x 2 V572T C ° S a+ 715 Sm " 



(75) 



,t n2 , / 3T 2 2493 . 2 

< (Jx) > = COS a +^f Sm 1 



(76) 



We can see that if a has been determined prior to the moment when the random walk in frequency begins to dominate, and its 
numerical value was not close to vr/2 or 3tt/2, then the information about this value is lost after (28607ra; 2 )/ (3/14) ~ 2995a; 2 //i4 
orbital revolutions. In the case when a is close to tt/2 or 3tt/2 the information will be lost about parameter x. Arranging 
observations in the interval before random walk begins a dominant source of the noise we can preserve information either 
about numerical value of a or x parameters. If the initial orbital phase becomes the non-informative parameter one has to 
average equations (p3|)-(|50|) with respect to a under assumption that distribution of probability for this parameter is uniform 
and concentrated in the interval from to 2tt (Bard 1974). 



7.7 Flicker Noise in Frequency Derivative (1// noise) 

The timing residuals are proportional to the forth power of the number of orbital revolution 



(77) 



The non-stationary part of the autocovariance function for flicker noise in frequency derivative can be expressed through the 
basic functions as (Kopeikin 1997b) 



16 



(78) 
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(85) 



Using (Eq) the covariance matrix for the non- stationary part of the flicker noise in frequency derivative is reduced to the form 



rr . 10 . . 10 . . 17 . . 17 . . 56 . . 

-Od al d b l - —d al d b 2 — —0 b id a 2 - O a ld b3 — ^p0blda3 + ^20 a 20 b 2 



(86) 
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Table 15. Elements of the covariance matrix M b of pulsar's parameters for flicker noise in frequency 
derivative. Quantities T = ttN and Q = cos T. The magnitude of noise /15 is omitted. 
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Table 16. Elements of the covariance matrix M ab of pulsar's parameters for flicker noise in frequency 
derivative. Quantities T = ttN and Q = cos T. The magnitude of noise h§ is omitted. 
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It is clear from (^) that the only elements of the matrix M^ b being not equal to zero are M^ b = M bl , M 2b = M b2 , and 
M^ b = M b3 , where the index 6=1, 14. In particular, the elements M^, M 22 , and M^ 3 are given 



ln(2T) 



613 
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(94) 



M± = h 5 T 2 



ln(2T) 
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(95) 



Mi 



ln(2T) + 
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(96) 



The matrix M ab is given in Tables (|l^)-(P^|). From these tables, equations (^)-(^), and discussions of previous sections we 
conclude that just the nicker noise in frequency becomes a dominant source of noise in timing residuals - the initial rotational 
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phase No, rotational frequency v, the initial orbital phase a, the semi-major axis x - are getting the non- informative parameters 
after a certain period of time which can be calculated in the same way as it has been done in previous sections. Determination 
of numerical values of parameters - u, x, and the orbital frequency n b can not be improved comparatively to those values 
which was obtained at preceding epoch when timing noise was not so red. This has a consequence for stability of the, so called, 
binary pulsar time (BPT) scale (Ilyasov, Kopeikin & Rodin 1998) which is discussed in section 8. Flicker noise in frequency 
may be produced by the stochastic gravitational wave background. Detailed discussion of the analysis of timing data from 
binary pulsars in the presence of this noise is given by Kopeikin (1997a). 

7.8 Random Walk in Frequency Derivative (l// 6 noise) 

The timing residuals are proportional to the fifth power of the number of orbital revolutions 

< ^ >= 765765^' < 97 > 
where is a constant quantity characterizing intensity of the noise. The non-stationary part of the autocovariance function 
for PN random walk is represented as (Kopeikin 1997b) 

R + (u,tj) = 1 Lh fi T 5 lMu)Mv) + ^iMu)Mv) + M")Mv)]+ (98) 



+^f2 WiW^W + ^W^W+^H^W] (99) 



5 5 

+ 7p3 h/>3(li)^2(«) + 1p2(u)lp 3 (v)} + T^1p3(u)lp 3 (v) (100) 



+ i2T^ hM u )^( w ) - 5 Mu)ip 5 (v) + wMu)Mv) (ioi) 



+10V'4(w)V , 3(«)-5V'5(w)V'2(v)+V6(w)V'l(«)]| ■ (102) 

Using formula (|2^) the covariance matrix for non-stationary part of PN is reduced to the form 

1 , ^5 



M - = To hsT 



5 5 20 

5 a i5 bl + — 8 al S b2 + —5 bl S a2 + -^^5 a 25 b2 (103) 



5 5 5 5 5 

+ 7^SalS b -s + -^^S bl Sa3 + <5 a 2 #63 + ^<5b2<5a3 + 7j^;f>a3<5b3 (104) 



^ 12T5 (^ a1 ^ 66 ^ bl ^ a6 ~~ 58 a 2&bh — 5&b2&ab + W8a3S b 4 + 10<5f,3<5 a 4) 



(105) 



The stationary part of the covariance matrix M~ b is given in Tables (p"7j)-(|l8|). Random walk in frequency is so strong at 
low frequencies that time derivatives of rotational phase also becomes the non-informative parameter after a certain period of 
time. Neither the orbital phase, nor the semi-major axis can be measured precisely. Variances of all parameters which depend 
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Table 17. Elements of the covariance matrix M , of pulsar's parameters for random walk in frequency 
derivative. Quantities T = ttN and Q = cos T. The magnitude of noise /ig is omitted. 
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Table 18. Elements of the covariance matrix M ab of pulsar's parameters for random walk in frequency 
derivative. Quantities T = ttN and Q = cos T. The magnitude of noise h§ is omitted. 
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on the numerical value of the initial orbital phase a should be averaged with respect to a assuming the uniform distribution 
of probability density function in the interval from to 2vr. 

8 KINEMATIC AND DYNAMIC PULSAR TIME SCALES 

In this section we discuss stability of time scales based on kinematic rotation of pulsars (PT scale) and orbital motion of 
pulsar in a binary system (BPT scale). This comprehensive analysis can be done due to using results of calculation of timing 
residuals and covaraince matrices given in the preceding sections. 

The methodology of applying pulsar timing data to fundamental metrology and time keeping service was explicitly 
formulated by Russian astronomers in 1979 (Shabanova et al. 1979, see also Ilyasov el al. 1989). Subsequent timing observations 
definitely proved consistency of this approach and made it clear that PT scale has a stability comparable with that of atomic 
clocks being placed in metrological centres (Ilyasov et al. 1989, Ginot and Petit 1991, Kaspi et al. 1994). However, it does not 
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mean that one should stop looking for other natural periodic physical phenomena which might be used for generating time 
scales being even more stable than PT. A new possibility to explore the problem is creation and maintenance of modified 
ephemeris time (ET) of classical astronomy based not on the orbital motion of planets in the Solar system (Guinot 1989) but 
on the orbital motion of pulsar in a close binary system (Petit and Tavella 1996, Ilyasov et al. 1998). It has been theoretically 
justified (Damour 1987, Kopeikin 1985, Schafer 1985) and observationally confirmed with accuracy 0.4% (Damour and Taylor 
1991, 1992) that orbital motion of pulsar in a binary system is governed by laws of General Relativity. For this reason one can 
predict the orbital motion of pulsar with extremely high precision on very long time intervals larger than 10 years generating 
in this way the new time scale BPT. However, the principal question arises about the limit on the stability of such time scale 
in the presence of red stochastic noise process in pulsar timing residuals. We answer this question in the present section. 

8.1 Binary pulsars as time keeping standards 

Standards of time and frequency have three important characteristics: (1) frequency stability, (2) reproduction of the time 
scale unit, and (3) the span of life time. Regarding these characteristics we note that the life time of binary pulsars reaches 
10 6 years and more which makes them very long-lived time standards. Orbital motion of pulsar in a close binary system is 
practically not subjected to external gravitational perturbations. For this reason the orbital motion has a very high intrinsic 
stability allowing reproduction of duration of one orbital revolution to very high precision. Binary pulsar time scale - BPT 
(Petit & Tavella 1996, Ilyasov et al. 1998) can be constructed according to the formula 

T = % + P b N + ^P b N 2 , (106) 

where To is an initial epoch for orbital phase, N is the number of orbital revolutions, P;, and Pi, are the orbital period and its 
time derivative reffered to the initial epoch To. This equation does not contain terms depending on high-order time derivatives 
of the orbital period as they equal identically to zero in the barycentric reference frame of the binary system (Kopeikin 1985, 
Kopeikin and Potapov 1994). However, radial acceleration of the binary system with respect to observer and its proper motion 
in the sky bring about appearance of the second and next high-order time derivatives both of Pb (Damour & Taylor 1992, Bell 
& Bailes 1996) and of the projected semi-major axis of the pulsar's orbit (Kopeikin 1996, 1997a) (see equation (Q) for more 
detail). These effects change the intrinsic values of Pi, and P;, and make more difficult their determination from observations. 

An additional serious problem arises with the possible presence in TOA the low-frequency noise of astrophysical origin. 
Indeed, intrinsic rotational motion of pulsar may be perturbed by the random walk in phase and/or its time derivatives 
(Cordes & Greenstein 1981), interstellar medium and relic stochastic gravitational waves change randomly propagation of 
electromagnetic pulses from the pulsar to observer being the reason for contamination of timing residuals by flicker noise 
(Mashhoon & Grishchuk 1980, Bertotti et al. 1983, Blandford et al. 1984, Backer & Foster 1990, Kopeikin 1997a). As a result, 
timing residuals represent a pure deterministic process upon which the stohastic additive noise is imposed. This noise consists 
of mixture of white noise of errors of observations and the low-frequency red noise of astrophysical origin. The spectrum S(f) 
of the red noise is described in Table [| It contains both stationary and non-stationary components (Kopeikin 1997b). The 
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presence of the red noise makes specific restrictions on the stability of both PT and BPT time scales which we briefly discuss 
in what follows. More detail on the subject may be found in the paper by Ilyasov, Kopeikin & Rodin (1998). 



8.2 BPT scale and its stability 

A complete treatment of the problem under consideration requires rather combersome calculations. In order to avoid mathe- 
matical difficulties we restrict ourselves to a simplified model of a binary system with the pulsar on the circular orbit so that 
all results of previous sections can be applied without any restrictions. Our goal is to get the optimal estimate of measured 
stability of BPT on the background of additive, low-frequency red noise in TOA. This problem had been partially considered 
(see, for instance, Van Trees 1968, Tikhonov 1983). Unfortunately, the methods developed in most previous publications can 
not be applied directly in our study so that we need to rely upon our own approach. In the problem under consideration 
the signal represents a linear combination of rotational phase of the pulsar, which is given by a polynomial of time, and 
quasi-periodic sinusoidal function with the argument being the orbital phase which is also a polynomial of time. As we have 
assumed the noise is additive to the signal, we can use the functional expression for the analyzed signal in the form: 

£(t) = M{t, /3a) + v e(t), 0<t<r, (107) 

where Af(t,/3 a ) is function of time describing deterministic part of the signal, depending on measured parameters (3 a , and 
given by equation (^); v is the pulsar rotational frequency; e(t) is a random process defined by equation (|B|) which differentiate 
between noises which are true instabilities of the pulsar clock, and apparent instabilities that arise from random galactic 
accelerations, gravitational waves, etc. The autocorrelation function of the noise process e(t) is given in Table ^, r is the total 
span of observational time. We also assume that observations are evenly spaced along the orbit with frequency m. It has 
allowed us to replace in all sums by integrals and to treat the stochastic process e(t) as continuous. 

For an adequate treatment of stability of time scales PT and BPT it is convenient to introduce two new parameters 
y — Su/u and v = Snt/nt, where 8v and 8n^ represent differences between true physical value of the parameter in question 
and its estimate obtained by the least squares method when fitting TOA of the pulsar. Variances of estimates of rotational 
and orbital frequencies (<jy, and, crjj, respectively) are convenient for comparing stability of two time scales - dynamic BPT 
and kinematic PT. Mathematical expressions for these variances are called true variances (Rutman 1978) and represent the 
close analogue of the Allan variance being used in fundamental metrology for analysing stability of time-frequency standards. 
Expressions for a v and cr y crucially depend on the type of noise present and can be evaluated from the variance of jii and from 
equation ( fli] ) and variances Pr, /3s. Their explicit dependence on the amplitude of noise h s , s — 0, 1, 2, . . . , 6 and duration of 
the total span of observations r are given in Table |l9| and are ploted schematically in Fig. ^ where the horizontal axis shows 
the total span of observations r and vertical axis numerates values of the parameters a y and a v . It is worth emphasizing 
that a y can not depend on orbital parameters of the pulsar because it characterizes instability of the intrinsic rotational 
frequency of the pulsar and, for this reason, can not be related with its orbital motion. On the other hand, the quantity a v 
naturally depends on the orbital period of the binary Pt = 2ir/nb and its projected semi-major axis x since these parameters 
characterize frequency and amplitude of the orbital sine wave in timing residuals. Additional functional dependence of a v 
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Figure 2. Characteristic behavior of variances characterizing stability of intrinsic rotational frequency of pulsar cr y and that of the 
orbital frequency of the binary system a v . In the case of evenly spaced observations theoretically predicted ratio a y /a v is equal to 
1.75tt x/P b ~ 1CT 3 -r 1CT 4 until white noise dominates in timing residuals. Displacement of two curves in subsequent moments of time 
depends on the relative magnitude of white and red noise components. 

on the initial orbital phase a in the case of random walk in phase is explained by the fact that the orbital frequency of the 
binary is one of the quadrature component of the sine function (the second quadrature component is the first time derivative 
of the projected semi-major axis x as follows from the definition of parameters fig and /3io given in Table 1). Low- frequency 
red noise with spectral index s = 2 distroys an equilibrium in the accuracy of simultaneous determination of the quadrature 
components and leads to the dependence of a v on the parameter a. In the case of red noise with index s > 3 the initial orbital 
phase is the non-informative parameter. Thus, we have averaged the variances with respect to a which, for this reason, does 
not appear in the corresponding expressions for er„. 
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In constructing Fig. || we have assumed that during the time interval [to,ti] white phase noise dominates. This noise is 
caused by the presence of errors of measurements in TOA and, possibly, by the uncorrelated noise in the intrinsic rotation of 
the pulsar. Various components of the low-frequency red noise commence to appear on time intervals r > n. Their amplitudes 
h s , as the experience of working with atomic clocks shows, gradually decrease as the spectral index s of the noise increases. 
For instance, in the model being considered in Fig. ^ the red nicker noise with spectral density S(f) ~ hi/ f dominates in 
timing residuals from the moment T\ till ti- Then, starting from the moment T2 up to T3 the random walk in phase with the 
spectral density S(f) ~ h^j '/ 2 dominates, and so on. In general, we assume that the longer time interval of observations the 
more significant the contribution of the red noise with higher numerical value of the spectral index s. Such behavior relates 
to the fact that usually the red noise with higher spectral index has lower intensity (that is, smaller amplitude h 3 ) and can 
give rise to a noticable contribution only on longer time interval r from the beginning of observation. 

It is interesting to note that behavior of (J v {t) is reminiscent of the, so-called, narrow-band frequency dispersion used in 
the fundamental metrology for characterizing stability of time-frequency standards when low-frequency noise is dominating 
(Rutman, 1978). However, in the timing model under consideration, a y (r), turns out to be significantly dependent on the 
non- stationary part of red noise with the spectral indices s = 3 or high (see Table As there is no satisfactory theoretical 
model of the non-stationary component such dependence of <t v {t) is an indication of necessity of working out the other, more 
practical approach for estimation of the variance characterizing stability of intrinsic rotational motion of pulsar and PT scale. 
It seems that the modified measure of the variance, which is called a z instead of a y , introduced by Taylor (1992) and based on 
the using of mathematical technique of transfer filter functions in spectral frequency domain is rather constructive and fruitful 
step ahead toward this direction. In particular, Matsakis et al. (1997) have recently derived the explicit expression for cr z (r) 
which is proportional to the second time derivative of pulsar's rotational frequency i>. We have used our analytic estimation 
of pulsar parameters in timing model (Q) in order to calculate u z . Results of the calculations is shown in Table [[9] along 
with a y and a v . We confirm the statement of Matsakis et al. (1997) that function <r z is independent from the non-stationary 
component of red noise. Moreover, it does not contain slow logarithmic trends explicitly appearing in a y for flicker noises with 
indices s = 3 and s = 5. 

Comparing the a z statistic with a v we find that a v is not sensitive to the non-stationary component of red noise as well 
as <t 2 , and does not contain the logarifmic trends. However, a v may depend on the initial orbital phase for red noises with 
indices s > 2 if span of observational time is not long enough to make that parameter to be non-informative. Our calculations 
also clearly show that measured value of the variance <t„(t) is not sensitive to red noise with index s < 2 and, for this reason, 
does not allow one to distingiush such red noise from the white one. Nevertheless, this specific behavior of a v (r) opens a 
possibility of using orbital motion of the pulsar as a new time reference being more stable than PT scale on longer intervals. 
Indeed, as one can see in Fig. ^, function ^(r)^ describing stability of PT scale begins to grow from the moment T3 but <j v (t) 
characterizing stability of BPT scale continues to decrease until the red noise with spectral index s > 5 begins to dominate. 

** Behavior of <r z is more representative for characterizing stability of PT scale. However, graphic behavior of a v looks approximately 
the same as that of a z . 
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Table 19. Dependence of variances <j v (t), <t z (t), and ct v (t) on the total span of observations r for red noise 
with different spectra. Herein, h s , (s = 1, 2 . . . , 6) is the amplitude of the corresponding noise component; 
Pt, is the orbital period and x is the projected semi-major axis of the binary orbit, a is the initial orbital 
phase, At is the interval of time between two successive observations, C:j,...,Ca are positive constants of 
order 0.05 - 0.5 reflecting dependence of the variances on the non-stationary component of the noise model. 
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This relative behavior of two variances is quite general and does not depend on the specific numerical value of intensities 
of red noise. Theoretical analysis reveals, that the minimum of a v (r) is reached significantly later than that of <J y (r) (and 
ct 2 (t) respectively). Depth of the minimum for cr y (T) (and cr 2 (r) depends on the noise with spectral index s = 3 which is 
produced by the large-scale inhomogeneities of the interstellar medium (Blandford et al. 1984). Depth of the minimum for 
u v (t) depends on the amplitude of noises with spectral indecies s = 5 which exist due to the stohastic background of relic 
gravitational radiation produced by physical processes in early universe (Mashhoon & Grishchuk 1980, Bertotti et al. 1983, 
Kopeikin 1997a). In principle, noise with spectral index s = 5 may arise also as a result of random fluctuations of the first 
time derivative of intrinsic rotational frequency of the pulsar (Cordes & Greenstein 1981). The existence of such fluctuations 
is extremely unlikely and, for this reason, they are not discussed in the following discussion. As an example we specify the 
numerical value of the minimum of ct z (t) for a single millisecond pulsar PSR B1937+21 which reach 10 -14 on time interval 2-3 
years. Existence of the minimum with subsequent turning-up of the curve a z (r) for this pulsar is explained by the dominance 
of random instabilities in the rotational phase of the pulsar (Kaspi et al. 1994). 

It is remarkable from principal and practical points of view that depth of the absolute minimum for the function, <j v (t), can 
be predicted with a high degree of certainty. It is defined by the amplitude of the relic stochastic gravitational wave background 
radiation with spectral index s = 5. More specifically, for binary pulsars with circular orbits the depth of minimum of a v (r) 
is defined by the expression: 

a v ~ 1.2 ■ W'™ y/n^PtX~\ (108) 

where Q g is energy density of stochastic gravitational wave background radiation per logarithmic unit interval of frequency, 
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Pi, is the orbital period of the binary measured in seconds, x is the projected semi-major axis of the pulsar's orbit measured in 



seconds, h is the Hubble constant in units 100 km c _1 Mpc _1 . It is worthwhile to point out (as the formula (108) shows) that 
numerical value of the minimum for a v (r) depends not only on the magnitude of fundamental cosmological parameters Q g and 
h but on the numerical values of orbital parameters of the binary system as well. Of course, for binary pulsars with elliptical 
orbits formula (5) will also include dependence on the eccentricity of the orbit. We shall consider this more complicated case 
elsewhere. 



Formula (108) yields the absolute (fundamental) value of the minimum for a v (r). However, the real minimum of a v (r) 
depends significantly on the amplitude of noise with spectral index s = 6. If the amplitude of this noise is less than that 
of the stochastic noise of gravitational wave backround radiation, then the real minimum of g v (t) coincides with absolute 



minimum given by formula (108). In the opposite case the minimum of a v (r) will be located above the absolute minimum 



( 105 ) . In any case, the observed numerical value of minimum for a v (r) may, in principle, reach the level 10~ J -i- 10~ under 
assumption that Q g h 2 < 10 -8 (Starobinsky 1979, Vilenkin 1981, Rubakov et al. 1982). The result obtained can be used in 
two ways. First, existence of fundamental minimum for <j v (t) allows to search for (or to establish the upper limit on) the 
stochastic gravitational wave radiation by means of using binary pulsars with long orbital periods and large ratio of Pb/x 
which is necessary for reducing the absolute value of minimum of <t v (t) as much as possible. That will reduce the interval of 
time r being required in order to reach the minimum. Second, using binary pulsars with short orbital periods and small ratio 
of Pb/x one can make the minimum of a v (j) as deep as possible in order to make the interval of stability of BPT scale as 
long as possible. We argue that determination of numerical value of minimum of a v (r) from observations may give very likely 
more precise indicator of the upper limit on the energy density of stochastic gravitational radiation than measuring the slope 
of function <j z (t). The matter is that determination of minimum of o"„(t) is more statistically reliable than one of the slope 
of <j z (t) which depends to a large extent on uncertainty in calculation of errors of measurement of the curve on long time 
intervals (Kaspi et al. 1994). 



9 CONCLUSIONS 

We have done analytic calculations of time dependence of residuals and covariance matrix of fitting parameters of a binary 
pulsar under assumption that observations are evenly spaced. The results have been used to estimate the observed stability of 
intrinsic rotational and orbital frequencies of the pulsar. It was shown that the dynamic time scale BPT may be more stable 
under certain circumstances than the kinematic PT scale. It was clarified that stability of BPT is restricted (like that of PT) 
by the presence of red noise even when the orbital motion of the pulsar can be considered as pure deterministic. Nevertheless, 
this restriction on stability of BPT is important only if observations are run on time intervals being significantly longer than 
characteristic time of instability of PT. This remarkable property of BPT leads us to hope that its usage may help in future 
to undertake deeper analysis of laws of gravitational physics. There is no doubt that long-term stability of BPT time scale 
makes it as a useful practical tool for fundamental metrology of time. 

Investigation of problem of construction and maintanace of BPT scale for studing various problems of modern astronomy 
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relates tightly to the problem of identification of nature and spectral characteristic of timing noise presenting in TOA. There 
is a hope that publication of papers of Deshpande et al. (1996) and D'Alessandro et al. (1997) devoted to detailed discussion 
of timing noise in the case of a sample of 18 single pulsars will attract a new interest to this problem. In particular, a reader 
of these papers will find a handful cases where the slopes are possibly more consistent with the spectral indices n — 3,5 etc. 
than say n = 4, 6. Of course, the uncertainty in some of these cases is large but still there is an indicative support for the 
flicker noise cases. 

Theoretical development of ideas relating to the study of stability of BPT time scale can be continued in several directions. 
Probably, one should consider binary pulsars with elliptical orbits having maximal ratio Pb/x which can improve quality of 
the time scale and, consequently, testing General Relativity in the strong field regime. In this problem it will be necessary 
to take into account possible variations of orbital elements of the binary caused by different, not well predicted factors, like 
flybys of stars or gravitational waves with frequncies being close to the orbital one (Chicone et al. 1998). It is not so difficult 
to show that for the BPT it is preferable to take binary pulsars with small enough ratio m p /m c (m v and m c are masses of 
pulsar and its companion respectively), long semi-major axis of the pulsar's orbit a p , short orbital period Pt, and sini ~ 1, 
that is large x = a p sin i/c. It would be desirable and worthwhile to create a network of reference binary pulsars and start 
their continuous timing observations in leading radio astronomical observatories of the world in order to make the concept of 
BPT time scale practically available as soon as possible. 
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